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I study numerically the problem of delamination of a thin film elastically attached to a rigid 
substrate. A nominally flat elastic thin film is modeled using a two-dimensional triangular mesh. 
Both compression and bending rigidities are included to simulate compression and bending of the 
film. The film can buckle (i.e., abandon its flat configuration) when enough compressive strain is 
applied. The possible buckled configurations of a piece of film with stripe geometry are investigated 
as a function of the compressive strain, ft is found that the stable configuration depends strongly 
on the applied strain and the Poisson ratio of the film. Next, the film is considered to be attached 
to a rigid substrate by springs that can break when the detaching force exceeds a threshold value, 
producing the partial delamination of the film. Delamination is induced by a mismatch of the relaxed 
configurations of film and substrate. The morphology of the delaminated film can be followed and 
compared with available experimental results as a function of model parameters. 'Telephone-cord', 
polygonal, and 'brain-like' patterns qualitatively similar to experimentally observed configurations 
are obtained in different parameter regions. The main control parameters that select the different 
patterns are the mismatch between film and substrate and the degree of in-plane relaxation within 
the unbuckled regions. 

PACS numbers: 



I. INTRODUCTION 

The use of a variety of coatings to enhance the perfor- 
mance of materials is widespread in many areas of science 
and technology. The deposition of the coating usually 
occurs at conditions (such as temperature, humidity of 
the environment, etc) very different of those to be found 
under work conditions. Due to this fact and to the dif- 
ferent nature of substrate and coating, typically large 
mismatch stresses appear between the film and the sub- 
strate. When the stresses are tensile on the film, this may 
fail duo the the nucleation of cracks that split the film to 
reduce to total mechanical energy of the system. A famil- 
iar example of this phenomenon is the cracking observed 
sometimes in paints and in the surface of mudi When 
mismatch stresses are compressive within the film, the 
most common mode of failure is called delaminationi2ii4 
the film partially detaches from the substrate to relieve 
the accumulated stress. Experimentally, different mor- 
phologies of the delaminated regions have been observed. 
In many cases delaminated regions of characteristic un- 
dulated geometries (referred to as 'telephone cords') ap- 
pear. 

Delamination can be understood on the basis of 
the well known phenomenon of buckling of elastic 
structures^ The first quantitative description of the 
buckling of a one dimensional elastic rod goes back to 
Euler. Upon compression, an elastic rod can reach a 
state in which the straight configuration is no longer the 
minimum energy configuration. At that point the rod 
acquires a non-straight configuration, namely, it buckles. 
For two dimensional elastic membranes the same kind of 
instability exists, but the problem becomes mathemati- 
cally much more involved. The equilibrium equations of 
a membrane were obtained by Foppl and von Karman 
(FvK)pi The FvK equations are two non- linear, coupled 



differential equations for the separation of a given point 
with respect to the flat configuration, and for the Airy 
potential for the in-plane displacement. The in-plane dis- 
placement can be reconstructed from auxiliary (linear) 
differential equations once the Airy potential is known. 

Delamination of a film that is elastically attached to a 
rigid substrate involves the interplay of two main ingre- 
dients. For a given initial form of a delaminated region 
(in which the interaction with the substrate is assumed 
to be absent) the FvK equations determine the geometry 
of the buckled configuration. The buckled film produces 
forces on the border that tend to detach the film further. 
In this way the film can continue its detaching by this 
mechanism, with the configuration of the buckled film 
adapting to the instantaneous form of the buckled re- 
gion. Thus, in general, buckling and delamination of the 
film must be solved self-consistently, and this is possible 
only numerically. 

A direct numerical solution of the FvK equations is 
rather difficult, because of its non-linear nature. This 
is why some people have considered the simulation of a 
more-or-less realistic film, consisting of atoms interact- 
ing via inter-atomic forces. The problem with this ap- 
proach is that in order to be able to model the bending 
of the film, a few layers of atoms are necessary, making 
the simulations rather inefficient, and suffering of impor- 
tant spurious effects associated to the anisotropy of the 
numerical lattice. 7 Here, to simulate the film I use the 
following strategy, that to a large extent eliminates the 
previously mentioned drawbacks. 8 The film is modeled on 
a physical basis by a triangular lattice formed by springs 
that account for the compression/stretching behavior of 
the material. In addition, bending springs are introduced 
in the sites of the triangular lattice that account for the 
bending rigidity of the film. The model can be shown 
to be elastically isotropic (in linear approximation), de- 
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scribing a material with Poisson ratio v = 1/3. Addi- 
tional bending springs are added to change the value of 
v in a controlled manner, keeping the system isotropic. 



By referring to the sketch in Fig. [Ja), each node % of 
the triangular lattice contributes with an elastic energy 
Hi of the form: 



where uy is the vector joining sites i and j, the index 
j labels the six neighbors of site i, in sequential order, 
j' in the second term is j + 3 (i.e, the neighbor opposite 
to neighbor j), and j" in the last term is j" = j + 1 
for j from 1 to 5, and j' = 1 for j = 6. The three 
terms of this expression are respectively the stretching 
energy, the bending energy, and the corrective term that 
is introduced if a film with v j= 1/3 is required. They are 
characterized by three different elastic constants, namely 
k s , kf,, and k u . The equilibrium configuration of a system 
described by the energy is a flat triangular lattice with 
lattice parameter Iq. 
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FIG. 1: (a) Detail of the triangular lattice that simulates 
the film, where the numbering of neighbors of site i and the 
equilibrium distance lo are highlighted, (b) A sketch of part of 
a partially delaminated film. Open symbols are the anchoring 
points in the substrate. The film is attached to the substrate 
to the left, and is detached to the right. Note the horizontal 
relaxation in the attached part, indicated by the thin lines, 
along which the spring constant kh is acting. 

From this atomistic expression for the energy of the 
system, the relation with averaged macroscopic quanti- 
ties can be evaluated. The result for the bending rigidity 
D, the two-dimensional compressibility B and the Pois- 
son ratio v of the bulk material of the film (which is 
assumed to have a thickness d) is& 
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The attaching (and eventual detaching) of the film to 
the substrate is the second crucial ingredient to include 



in order to obtain a realistic behavior of the system. One 
possibility (as used for instance in Rcf. 7 ) is to consider 
that when attached to the substrate the film is rigidly 
pinned to it (i.e., it cannot move at all) whereas, if the 
force exerted by the substrate onto the film exceeds some 
threshold value, the interaction to the substrate is set 
irreversible to zero. Experimentally, it has been em- 
phasized that detaching by perpendicular separation has 
usually a much lower threshold than detaching by shear. 4 
This suggests that an appropriate situation is to consider 
that detach occurs when the perpendicular force between 
film and substrate exceeds some maximum value f m ax- 
This was the scheme adopted ml. We will see that this 
prescription does not suffice to reproduce the large phe- 
nomenology that is observed experimentally. 1 think that 
the main lacking ingredient is the possibility of in-plane 
relaxation of the film in the case it is still attached to 
the substrate. This relaxation is introduced here in the 
energetics of the model in the following way (see Fig. 
Elb)). Let us consider the substrate to be located at 
z = 0. For each point of the film there is an anchor- 
ing point on the substrate to which this is attached to. 
The coordinates of the anchoring point are (xq, ya, 0). 
The coordinates of the corresponding point in the film 
are (u x , u y , 0). The contribution to the energy from this 
link is \kh[{u x — xq) 2 + (u y — yo) 2 }- The value of kh deter- 
mines the extent of in-plane relaxation. For kh — > oo we 
recover the rigid attaching condition. When the z force 
that is necessary to apply to a given point in the film 
to keep it at zero z coordinate overpasses the threshold 
value fmax, the value of the corresponding kh is set ir- 
reversibly to zero, and the z coordinate is not any more 
restricted to be zero. The anchoring points in the sub- 
strate form a triangular lattice, corresponding to that of 
the film. However the lattice parameter of this lattice Iq 
is assumed to be smaller than the one Iq in the film. This 
is the way in which a mismatch between film and sub- 
strate appears, and the origin of the tendency of the film 
to delaminate. Macroscopically, the attaching to the sub- 
strate can be quantified by an elastic coefficient we call 
E s , that is defined as the ratio between the applied hor- 
izontal force per unit area of the film and the horizontal 
displacement this force produces. The relation between 
kh and E s is easily found to be 
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In addition to the energetics, the time evolution of the 
system must be defined. This evolution will be assumed 
to be overdamped, namely it is assumed that kinetic en- 
ergy plays no role in the evolution. Then the positions 
of the nodes of the lattice evolve according to: 



du _ X SH 
dt <5u 



thus A sets the time scale for the evolution. 

It is also convenient to introduce non-dimensional pa- 
rameters to describe the simulations. From Q and JIJ, 
it can be seen that the simulations can be parametrized 
by the values of k b /(l 2 k s ), k u /(lok s ), and f ma x/(k s lo)- A 
non-dimensional time t can also be introduced through 
t = Xk s t. All simulations presented below where made in 
systems with periodic boundary conditions, using a first 
order method to integrate QJ, with a time step St = 0.3. 



II. BUCKLING PATTERNS OF A STRIPE 

Before studying the full interplay between buckling 
and delamination, it is convenient to consider the buck- 
ling patterns obtained for a fixed dclaminated region. I 
present the case of a delaminated region with stripe ge- 
ometry because this case has been considered previously 
in the literature^ and then it will be useful to validate 
the numerical model. In addition we will see that already 
in this simple case the behavior is non-trivial. Thus we 
consider a stripe piece of the film of width w, that is de- 
tached from the substrate. The rest of the film is assumed 
to be rigidly attached to the substrate (i.e., kh — ► oo). 
In this way, the only effect of this part of the film is to 
provide a clamped boundary condition for the stripe. In 
addition, due to the presence of the substrate, the stripe 
can buckle only toward positive values of z. Under these 
conditions we will look for the stable configuration of the 
stripe for different parameters. 

We define the mismatch strain between film and sub- 
strate e as e = (1 — Iq/Io)- When the stripe is loaded 
isotropically in compression (i.e, e > 0), the film remains 
flat until a critical value E\ is reached. At this point a uni- 
form wrinkle along the longitudinal direction is formed. 
This configuration is usually refer to as the Euler column. 
The value of E\ can be analytically calculated, and it is 
known to beiS 
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An equivalent expression can be obtained by expressing 
D and B in terms of the Young modulus of the material 
E, es indicated in 9 . The result is 
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where we see that the value of E\ is mainly controlled by 
the ratio between thickness and width of the delaminated 
stripe. 



The value of Ex can a ls° be expressed in terms of the 
parameters of the model by using (J2J , as 
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(4) where n is the number of rows of particles that compose 



the delaminated stripe (i.e., w 
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., mi"/, and where the 
previously defined non-dimensional parameter ap- 
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pears. 

The Euler column configuration mostly relieves the 
perpendicular stress on the stripe, but leaves an im- 
portant degree of longitudinal compression in it. If the 
isotropic compression e is increased, a secondary buckling 
is expected at some mismatch strain £2 , in which the lon- 
gitudinal uniformity of the Euler column is lost, and one 
of two different patterns is expected*^ in one case (occur- 
ring if the Poisson ratio v of the film is larger than a value 
of approximately 0.255) the Euler column transforms in 
an undulated wrinkle. The second case corresponds to 
a pattern (called sometimes 'varicose') appearing if v is 
lower than this value, in which the uniform maximum 
of the Euler column is distorted and transforms in a se- 
quence of small hills and shallow valleys. The varicose 
pattern preserves the symmetry of the configuration with 
respect to the axis of the delaminated region, whereas the 
undulated wrinkle does not. These results have been ob- 
tained analytically doing a perturbative analysis of the 
FvK equations considering small distortions of the Eu- 
ler column. The analysis shows alsoiS that £2/61 is a 
function of the Poisson ratio of the film only. Using the 
numerical model presented here, we can test this predic- 
tion, and in addition we can explore arbitrary distortions 
of the film. 

Numerical results in the stripe geometry were obtained 
taking one of the main directions of the triangular simula- 
tion lattice along the stripe. The stripe was described by 
24 rows of this triangular mesh, and the value of k b j(l\k s ) 
was taken to be 0.05iii Results are synthesized in the 
phase diagram of Fig. |2J and in Figs. and [3 where 
the obtained configurations are shown for different values 
of the two main control parameters e/e\ and v. In agree- 
ment with the analytical predictions, we can see that by 
increasing e, the instability of the Euler column corre- 
sponds to an non-symmetric pattern for v > v c ~ 0.25, 
whereas it corresponds to a symmetric one for v < v c . 
However, we also see that when larger values of e are 
considered, the stability range of the symmetric pattern 
is (at least for v > 0.05) only a rather small wedge in 
the phase diagram. For sufficiently large values of e the 
undulated wrinkle pattern is observed always to be the 
stable one. 

The symmetric and non-symmetric sectors of the phase 
diagram in Fig. [2] have been divided by a dotted and 
dashed line, respectively. I describe now their meaning. 
The dashed line in the non-symmetric sector marks only a 
crossover: At the right of this line, the undulated wrinkle 
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FIG. 2: The v-e parameter space, with the different sectors 
in which different qualitative behavior is observed. The lines 
separating different regimes were obtained by doing indepen- 
dent simulation at different point of the e-v plane. Some of 
them are presented in the next three figures. See the text 
and Figures and|^| Simulations were performed using a 
triangular lattice of 200 sites along the stripe, and 24 rows in 
the perpendicular direction. The dimensionless combination 
of parameters kb/(k s lo) was taken to be 0.05. The value of 
the critical strain ei is in this case ei ~ 0.0032. 



is well develop, in the sense that it has a maximum of al- 
most constant amplitude that undulates along the stripe, 
see for instance Fig. 0^e). At the left of the dashed lined, 
there is a clear modulation of the amplitude (Fig. EI C ))- 
I emphasize again that the transition between the two 
different behaviors is simply a crossover. 

The dotted line within the symmetric sector in Fig. 
121 marks an interesting behavior change in the region of 
small or negative values of v. To the left of the dotted 
line we have the symmetric pattern as described analyt- 
ically in RefiiS. To the right of the dashed line, this 
pattern gets modified by the appearance of deep valleys 
interrupting the normal sequence of small hills and shal- 
low valleys (see Fig. OJc)). Note that in the center of the 
deep valleys, the film is almost touching the substrate. 
The appearance of these deep valleys is abrupt, namely 
they do not grow deeper and deeper as some parameter 
is changed, but they simply exist or do not exist. The 
periodicity of the appearance of deep valleys is difficult 
to obtain from the simulations, but the information col- 
lected from different simulations suggests that close to 
the dotted line of the phase diagram deep valleys are 
well separated from each other, whereas when moving to 
larger values of e they become closer, eventually leaving 
no shallow valleys in between. A qualitative explanation 
for the appearance of these deep valleys could be based 
on the following observation. In a material with negative 
Poisson ratio, a compression in one direction is accom- 
panied by a compression in the perpendicular direction. 
Transforming a shallow valley of a varicose pattern into a 
deep valley produces a perpendicular compression of the 
film, that implies a tendency to compress also in the lon- 
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FIG. 3: (Color online) The buckled configuration of the stripe, 
for v — 1/3, and e/ei = 6.25, 7.5, 9.4, 12.5, and 25, respec- 
tively form (a) to (e). The gray (color) scale indicates the 
departure of the film from the substrate, in units of the dis- 
cretization distance Iq. 



gitudinal direction. In this way part of the longitudinal 
compression is relieved, and this can stabilize the deep 
valley. 

The previous results show that, for the usual cases of 
materials with positive Poisson ratio close to 1/3, the 
undulated wrinkle pattern is the most likely to be ob- 
servable in experimental situations in which stripe de- 
laminated regions occur. 



III. DELAMINATION 

In this Section I consider the full interplay of buckling 
and delamination, and present the different kinds of de- 
lamination morphologies that can be obtained with the 
present model. From now on I restrict to the case of a 
film with v = 1/3 (i.e., k u = 0). When buckling and 
delamination both occur, there are two natural length 
parameters in the problem whose meaning is convenient 
to emphasize from the beginning. The first length scale 
Li is the typical size of a blister that is in equilibrium 
upon further delamination. AN estimation of L\ goes 
as follows. For a buckled region of typical linear size L\ 
(L\ may be the for instance the diameter of a circular 
blister, or the width of an Euler column) under a defor- 
mation e, the typical curvature of the buckled film at 
the points of contact with the substrate is ~ y/e/Li^ 
The torque T per unit length exerted by the delaminated 
film on the part attached to the substrate is obtained by 
multiplying this curvature by the bending rigidity D, i.e, 
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FIG. 4: (Color online) Same as Fig. H for v = 0.075, and 
e/ei — 4.0, 4.7, 5.6, 6.25, and 6.9, respectively form (a) to 
(e). 
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FIG. 5: (Color online) Same as Fig. [3]for v = 0, and e/si 
3.12, 3.75, 4.37, 6.25, and 6.9, respectively from (a) to (e). 



T = \feD I ' L\. The equilibrium condition for the film will 
correspond to this torque been equal to a critical max- 
imum value that the interaction with the substrate can 
sustain. I note this maximum torque as T max . Thus we 
obtain the order of magnitude of L\ as L\ ~ \/sD/T max . 
It is important to emphasize that the condition of a max- 
imum torque per unit length that I use here is totally 
compatible with the condition of a maximum force used 
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FIG. 6: (Color online) Temporal evolution of the delaminated 
regions for e = 0.09, k b /(k s l%) = 0.05, / m<M /(Wo) = 0.03, 
and kh = k 3 . Note that the gray (color) scale goes from 
darker to brighter when the vertical coordinate of the film 
increases, but non-delaminated regions are plotted as white 
for clarity. Vertical scale is in units of lo- Dimensionless time 
t for the four views are r = 5 x 10 3 , 2 x 10 4 , 3.5 x 10 4 , and 
2 x 10 5 . 



to describe the model in the previous Section. In fact, 
note first of all that a torque per unit length has the cor- 
rect units of a force. In addition, a detailed consideration 
of the geometry of the triangular lattice used, shows that 
T m ax hi the continuum description and f m ax of the atom- 
istic description are simply related by T max — fmaxVi/2. 
Then, in the microscopic units of the model L\ can be 
written also as L\ ~ \fekbl ' fmax- 

A second characteristic distance is related to the exis- 
tence of horizontal relaxation within the unbuckled film. 
In fact, consider a film uniformly attached to the sub- 
strate. If we force a horizontal displacement of some 
point of the film, this perturbation will decay away in 
a typical distance L2, which depends on the ratio be- 
tween the stiffness of the film and the stren gth of the 
interaction to the substrate, namely L2 ~ loy/k s /kh, or 
in macroscopic units L 2 ~ y/Bd/E s . The ratio LijL\ 
is a measure of the importance of horizontal relaxation. 
This will be negligible if L2/L1 < 1, or important if 
Li/Li > 1. 

In order to observe the delamination process, I seed 
the simulation with a configuration in which a small part 
of the film is detached. This portion was chosen as a 
small rectangular piece, mis-oriented with respect to the 
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FIG. 7: (Color online) Same as Fig. ©for e = 0.085, and k h = 
k s , at times r = 5 x 10 3 , 2 x 10 4 , 3.5 x 10 4 , and 2 x 10 5 . Note 
in particular that the last panel shows a stable configuration 
in which the film does not delaminate further. 
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FIG. 9: (Color online) Same as Fig. ©for e = 0.07, k h = 
fcs/100, at times r = 2.4 x 10 3 , 5.4 x 10 3 , 3.5 x 10 4 , and 
1 x 10 5 . 
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numerical mesh to minimize spurious effects associated 
to it . In a first stage of the simulation I take f max — > oo. 
In this way buckling of the seeded delaminated region 
occurs, but the film does not delaminate further. After 
a stationary configuration is achieved, f ma x is put to a 
finite value, time is reset to zero, and the evolution of 
the film is calculated. In all cases presented the numer- 
ical lattice has a total of 300 x 300 nodes, and periodic 
boundary conditions are used. 

I present first some result in the case the in which hor- 
izontal relaxation is small, namely k^ = k s . This corre- 
sponds to a distance of the order of the discretization 
distance Zo , i.e., a very small value. In fig. [S]we see a se- 
quence of configurations corresponding to e = 0.09. We 
see how in this case the film delaminates almost com- 
pletely. This kind of pattern was called 'brain-like' in^. 
In fig. U| we see the situation for a slightly lower value of e, 
namely e = 0.085. After the delamination of a small re- 
gion near the original defect, the evolution stopped com- 
pletely, and delamination halted. We see that there is a 
rather sharp transition between a blister that is not able 
to grow for low strains, to an almost completely delami- 
nated film for larger strains. 

In the case in which the value of horizontal relaxation 
is higher, there are important differences in morphology. 
Figs. 00 E3 and 1 101 show the evolution of the same original 
defect for the case of kh = k s /100 (which corresponds to 
a relaxation distance L2 approximately ten times Zo) and 
e = 0.09, 0.07, and 0.065, respectively. For large values 
of the strain, the delaminated region propagates to all 
the film, however, contrary to the previous case, there 
are large sectors in which the film remains attached to 
the substrate. The typical size of these sectors is roughly 
determined by L2, whereas the typical width of the de- 
laminated stripes between sectors is order L\. Note that 
the propagation occurs via the formation and branching 
of undulated blisters, which at the end form a pattern 
of polygons bounded by undulated wrinkles. When e is 
reduced, the branching of the undulated blisters reduces, 
giving rise to a pattern like that in Fig. In a narrow 
band of values of the applied strain (Fig. ITUI) . a single 
one dimensional delaminated region (a 'telephone cord') 
may occur. For lower values of e (not shown) the original 
blister is not able to grow at all. 

IV. CONCLUSIONS 

I have presented results about the buckling of a thin 
film attached to a substrate. Buckling is induced by a 



mismatch in the relaxed configurations of film and sub- 
strate. For a fixed geometry of the delaminated region, 
namely a uniform stripe, the buckling patterns were ob- 
tained as a function of the strain mismatch e and the 
Poisson ratio of the film v. The numerical results confirm 
the findings of previous analytic treatments about the in- 
stability mode of the uniform Euler column: a undulated 
wrinkle for v > 0.25, and a symmetric, 'varicose' pattern 
for v < 0.25. However, the present numerical analysis, 
which is not restricted to small distortions of the Euler 
column, shows that the stability of the varicose pattern, 
in the region of experimental interest close to v ~ 1/3, is 
restricted to a narrow interval of e values, then limiting 
the possibility of observing this pattern experimentally. 

In the second part of the paper I considered the in- 
terplay of buckling and delamination, and argued that a 
key factor determining the morphology observed is the 
possibility or not of horizontal relaxation in the non- 
delaminated part of the film. For small horizontal relax- 
ation, and large enough e, the film delaminates almost 
completely, generating a 'brain-like' structure. When e 
is reduced, there is an abrupt transition to a situation in 
which the delaminated region is not able to grow. 

For the case of large horizontal relaxation in the non- 
delaminated part of the film, at large e the growing occurs 
via elongation and branching of undulated delaminated 
regions. At the end a polygonal pattern delimited by 
undulated stripes is observed. The main characteristic 
of this pattern is the existence of regions in which the 
film does not delaminate. When e is reduced, the ten- 
dency of the undulated stripes to branched is reduced, 
and eventually a regime where a single undulated stripe 
(a 'telephone-cord') grows is observed. If s is reduced 
further the delamination stops completely. 

The simulations presented clarify the origin of the dif- 
ferent morphologies observed experimentally. I want to 
emphasize in particular the importance of horizontal re- 
laxation in obtaining different kinds of patterns, most re- 
markably the necessity of horizontal relaxation to obtain 
the well known telephone-cord pattern. 
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